Julia

Paquetes para matemáticas

David Gómez-Castro

March 1, 2025

Entre otras ventajas, los algoritmos matemáticos en julia pueden quedar escritos de manera muy agradable

https://github.com/mossr/BeautifulAlgorithms.jl

LinearAlgebra.jl

Es parte de Standard Library

Operaciones básicas

using LinearAlgebra
A = [1 2 3; 4 1 6; 7 8 1]
3×3 Matrix{Int64}:
 1  2  3
 4  1  6
 7  8  1
@show tr(A)
@show det(A)
println("inv(A) = ")
inv(A)
tr(A) = 3
det(A) = 104.0
inv(A) = 
3×3 Matrix{Float64}:
 -0.451923   0.211538    0.0865385
  0.365385  -0.192308    0.0576923
  0.240385   0.0576923  -0.0673077

La sintaxis para sistemas es la habitual

b = [1,2,3]
A \ b 
3-element Vector{Float64}:
 0.2307692307692308
 0.15384615384615383
 0.15384615384615385

Y las descomposición espectral

eigvals(A) 
3-element Vector{Float64}:
 -6.214612641961068
 -1.554026596484783
 10.768639238445848

Matrices especiales

Type   Description
Symmetric Symmetric matrix
Hermitian Hermitian matrix
UpperTriangular Upper triangular matrix
UnitUpperTriangular Upper triangular matrix with unit diagonal
LowerTriangular Lower triangular matrix
UnitLowerTriangular Lower triangular matrix with unit diagonal
UpperHessenberg Upper Hessenberg matrix
Tridiagonal Tridiagonal matrix
SymTridiagonal Symmetric tridiagonal matrix
Bidiagonal Upper/lower bidiagonal matrix
Diagonal Diagonal matrix
UniformScaling Uniform scaling operator
A = UnitUpperTriangular( [ 1 2 3 ; 1 2 3 ; 1 2 3 ] )
3×3 UnitUpperTriangular{Int64, Matrix{Int64}}:
 1  2  3
 ⋅  1  3
 ⋅  ⋅  1

Más sobre matrices especiales

Factorizaciones

lu permite exportar L, U y vector de permutaciones, o matriz de permutaciones

A = rand(3,3)
B = lu(A)
LU{Float64, Matrix{Float64}, Vector{Int64}}
L factor:
3×3 Matrix{Float64}:
 1.0        0.0        0.0
 0.13597    1.0        0.0
 0.0384478  0.0634868  1.0
U factor:
3×3 Matrix{Float64}:
 0.672019  0.937098  0.477768
 0.0       0.551104  0.0154662
 0.0       0.0       0.127688

Julia incluso puede elegir una “buena” factorización para nuestro problema

factorize(A)
LU{Float64, Matrix{Float64}, Vector{Int64}}
L factor:
3×3 Matrix{Float64}:
 1.0        0.0        0.0
 0.13597    1.0        0.0
 0.0384478  0.0634868  1.0
U factor:
3×3 Matrix{Float64}:
 0.672019  0.937098  0.477768
 0.0       0.551104  0.0154662
 0.0       0.0       0.127688
factorize(A'*A)
Cholesky{Float64, Matrix{Float64}}
U factor:
3×3 UpperTriangular{Float64, Matrix{Float64}}:
 0.678694  1.02194   0.489495
  ⋅        0.547024  0.0228474
  ⋅         ⋅        0.127375

Volveremos sobre la importancia de factorizar en la próxima sesión.

Más sobre sistemas lineales

Para sistemas “difíciles” o cuando queramos la máxima velocidad, podemos utilizar el paquete LinearSolve.jl

using LinearSolve

A = rand(4, 4)
b = rand(4)
prob = LinearProblem(A, b)
sol = solve(prob)
sol.u
4-element Vector{Float64}:
 -0.16499821190175334
  0.5894600876191373
  0.7559934725184622
 -0.09049747220920887

Implementa muchos resolvedores optimizados y tiene, además, un sistema de selección de solvers adecuados.

Video explicativo en Youtube

Estadística y teoría de la medida

Random.jl es parte de Standard Library

using Random 
@show rand()
@show rand( (-1,2) )
@show rand(5:10);
rand() = 0.7548164892900898
rand((-1, 2)) = 2
rand(5:10) = 10

Statistics.jl es parte de Standard Library.

Es reseñable la familia JuliaStats: DataFrames, Distributions, HypothesisTests, TimeSeries, …

Distributions.jl

The Distributions package provides a large collection of probabilistic distributions and related functions. Particularly, Distributions implements:

  • Sampling from distributions
  • Moments (e.g mean, variance, skewness, and kurtosis), entropy, and other properties
  • Probability density/mass functions (pdf) and their logarithm (logpdf)
  • Moment-generating functions and characteristic functions
  • Maximum likelihood estimation
  • Distribution composition and derived distributions (Cartesian product of distributions, truncated distributions, censored distributions)
using Random, Distributions
d = Normal()
d + 1.0
Normal{Float64}(μ=1.0, σ=1.0)
@show mean(d)
@show rand(d, 3)
mean(d) = 0.0
rand(d, 3) = [0.02207721830462669, -0.6724539961700556, 2.254074922578876]
3-element Vector{Float64}:
  0.02207721830462669
 -0.6724539961700556
  2.254074922578876
x = [0, 0.4, -0.4] 
fit(Normal, x)
Normal{Float64}(μ=0.0, σ=0.32659863237109044)

Ver también MeasureTheory.jl

Turing.jl

Turing is a general-purpose probabilistic programming language for robust, efficient Bayesian inference and decision making. Current features include:

  • General-purpose probabilistic programming with an intuitive modelling interface;
  • Robust, efficient Hamiltonian Monte Carlo (HMC) sampling for differentiable posterior distributions;
  • Particle MCMC sampling for complex posterior distributions involving discrete variables and stochastic control flow; and
  • Compositional inference via Gibbs sampling that combines particle MCMC, HMC and random-walk MH (RWMH).

Un ejemplo sencillo

Plots.jl

Tutorial básico

Creando animaciones

using Plots 
anim = @animate for n=1:20
    plot( x-> sin(n*x) , label="")
end;

gif(anim, "animation.gif",  fps = 30)

mp4(anim, "video.mp4",  fps = 30)

En un notebook se puede hacer

using Plots
@gif for n=1:20
    plot( x-> sin(n*x) , label="" )
end

Pintando polígonos

x,y =  [0.0, 0.5, 1.0],  [0.0, 1.0, 0.0]
plot(x[[1:end;1]] , y[[1:end; 1]], fill = 0)

Otras librerías de representación

SciML

SciML is the combination of scientific computing techniques with machine learning

Ecuaciones Diferenciales

Resolución de EDOs

using DifferentialEquations
f(u,p,t) = 1.01*u
u0 = 1/2
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan)
sol = solve(prob, Tsit5(), reltol=1e-8, abstol=1e-8)

using Plots
plot(sol,linewidth=5,title="Solution to the linear ODE with a thick line",
     xaxis="Time (t)",yaxis="u(t) (in μm)",label="My Thick Line!") # legend=false
plot!(sol.t, t->0.5*exp(1.01t),lw=3,ls=:dash,label="True Solution!")

julia implementa por defecto Tsit5 en lugar de Dormund-Prince

ModelingToolkit.jl

ModelingToolkit.jl is a modeling language for high-performance symbolic-numeric computation in scientific computing and scientific machine learning. It then mixes ideas from symbolic computational algebra systems with causal and acausal equation-based modeling frameworks to give an extendable and parallel modeling system. It allows for users to give a high-level description of a model for symbolic preprocessing to analyze and enhance the model. Automatic symbolic transformations, such as index reduction of differential-algebraic equations, make it possible to solve equations that are impossible to solve with a purely numeric-based technique.

ModelingToolkit.jl una capa más sobre DifferentialEquations.jl

using ModelingToolkit, Plots
using ModelingToolkit: t_nounits as t, D_nounits as D

@mtkmodel FOL begin
    @parameters begin
        τ = 3.0 # parameters
    end
    @variables begin
        x(t) = 0.0 # dependent variables
    end
    @equations begin
        D(x) ~ (1 - x) / τ
    end
end

using OrdinaryDiffEq
@mtkbuild fol = FOL()
prob = ODEProblem(fol, [], (0.0, 10.0), [])
sol = solve(prob)

plot(sol)

Vídeo con un ejemplo industrial

JumpProcesses.jl

FeNICS.jl

FEniCS.jl is a wrapper for the FEniCS library for finite element discretizations of PDEs. This wrapper includes three parts:

Installation and direct access to FEniCS via a Conda installation. Alternatively one may use their current FEniCS installation. A low-level development API and provides some functionality to make directly dealing with the library a little bit easier, but still requires knowledge of FEniCS itself. Interfaces have been provided for the main functions and their attributes, and instructions to add further ones can be found here. A high-level API for usage with DifferentialEquations. An example can be seen solving the heat equation with high order adaptive timestepping.

Ecuaciones en Derivadas Parciales

Gridap.jl

Gridap provides a set of tools for the grid-based approximation of partial differential equations (PDEs) written in the Julia programming language. The main motivation behind the development of this library is to provide an easy-to-use framework for the development of complex PDE solvers in a dynamically typed style without sacrificing the performance of statically typed languages. The library currently supports linear and nonlinear PDE systems for scalar and vector fields, single and multi-field problems, conforming and nonconforming finite element discretizations, on structured and unstructured meshes of simplices and hexahedra.

It has a very reach library of worked examples

JuliaFEM

The JuliaFEM project develops open-source software for reliable, scalable, distributed Finite Element Method.

Optimización: JuMP.jl

JuMP is a domain-specific modeling language for mathematical optimization embedded in Julia. It currently supports a number of open-source and commercial solvers for a variety of problem classes, including linear, mixed-integer, second-order conic, semidefinite, and nonlinear programming.

JuMP is a package for Julia. From Julia, JuMP is installed by using the built-in package manager.

import Pkg
Pkg.add("JuMP")

You also need to include a Julia package which provides an appropriate solver. One such solver is HiGHS.Optimizer, which is provided by the HiGHS.jl package.

import Pkg
Pkg.add("HiGHS")

Ejemplo

Limitaciones y alternativas

Even if your problem is differentiable, if it is unconstrained there is limited benefit (and downsides in the form of more overhead) to using JuMP over tools which are only concerned with function minimization.

Otras opciones

Cálculo simbólico

Symbolics.jl

using Symbolics
@variables x y

\[ \begin{equation} \left[ \begin{array}{c} x \\ y \\ \end{array} \right] \end{equation} \]

z = x^2 + y

\[ \begin{equation} y + x^{2} \end{equation} \]

simplify(2x + 2y)

\[ \begin{equation} 2 \left( x + y \right) \end{equation} \]

Derivadas

@variables x g(x)
Dt = Differential(x)
Differential(x)

. . .

expression = expand_derivatives(Dt(x * g))

\[ \begin{equation} g\left( x \right) + x \frac{\mathrm{d} g\left( x \right)}{\mathrm{d}x} \end{equation} \]

. . .

expr2 = substitute(expression, Dict([g => x]))

\[ \begin{equation} x + x \frac{\mathrm{d}}{\mathrm{d}x} x \end{equation} \]

. . .

expand_derivatives(expr2)

\[ \begin{equation} 2 x \end{equation} \]

### Limitaciones

There is a list available of known missing features.

SymPy.jl

SymPy is a Python library for symbolic mathematics.

This package is a wrapper using PyCall

With the excellent PyCall package of julia, one has access to the many features of the SymPy library from within a Julia session.

It requires Python to run in the background. To install it sometimes it is needed to make a few

using Pkg
Pkg.add("Conda") #  if needed
using Conda
Conda.update()

The only point to this package is that using PyCall to access SymPy is somewhat cumbersome.

Symbolics.jl vs Sympy.jl

Redes neuronales

https://commons.wikimedia.org/wiki/File:Colored_neural_network.svg

Una capa

https://commons.wikimedia.org/wiki/File:Single_layer_ann.svg
  • \(x_1, x_2, \cdots\) son las entradas

  • \(y_1, y_2, \cdots\) son las salidas

  • \(w_{11}, w_{12}, \cdots\) son los

  • Se calcula \(\displaystyle b_1 = w_{11} x_1 + w_{21} x_2 + w_{31} x_3 + \cdots\)

  • Y luego \(y_1 = f ( b_1 )\)

  • \(f\) es la llamada

Denotamos \(\mathbf x = (x_1, \cdots, x_p)\), \(\mathbf y = (y_1, \cdots, y_q)\), \(\mathbf W = (w_{11}, w_{12}, \cdots )\)

Red de neuronas

Se pueden hacer redes

de varias capas:

la salida de una

es la entrada de la otra.

Si tenemos muchos pares de entradas y sus respectivas salidas

\[ \widehat{\mathbf y_i} = F(\mathbf x_i, \mathbf W) \]

Buscar los mejores pesos \(\mathbf W\) se llama

entrenar la red

Hay muchas generalizaciones de esta idea: redes profundas (muchas capas), redes convolucionales, …

¿Cómo de potente es esto? Vídeo de Google DeepMind

Evaluación de la red: forward-propagation

Al algoritmo para, dados los pesos y la entrada calcular la salida se lo llama forward propagation.

Entrenamiento de redes

Se separan los datos en un conjunto de entrenamiento y conjunto de pruebas

Para ajustar la red tenemos que elegir una función de pérdida \(d({\mathbf y}, \widehat{\mathbf y})\). Por ejemplo \(|{\mathbf y} - \widehat{\mathbf y}|^2\).

A partir de ella elegimos una función de coste, que tenga en cuenta todos las pérdidas. Por ejemplo el error cuadrático medio \[ J(W) = \frac{1}{N} \sum_{i=1}^N |{\mathbf y_i} - \widehat{\mathbf y_i}|^2 \]

Así:

  • Se elige \(\mathbf W\) minimizando el coste sobre el conjunto de entrenamiento
  • Se comprueba que es “bueno” en general el coste con este \(\mathbf W\) sobre el conjunto de prueba (que debe ser distinto al de entrenamiento)

Hay múltiples opciones para la elección de función de coste y hacer la optimización eficientemente es gran parte de la dificultad.

Optimización en el entrenamiento

Para optimizar se utilizan múltiples procedimientos.

El ejemplo más sencillo es el descenso por gradiente \[ W_{n+1} = W_n - \gamma \Delta J (W_n). \]

Para calcular \(\nabla J\) se emplea la regla de cadena. A este proceso lo llama back-progation.

Un ejemplo detallado en Julia

Flux.jl

Flux es una de las mejores librerías de ML en Julia.

Video de introducción

Ajustar una línea

Un ejemplo no trivial

Otros paquetes

MLJ.jl